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Abstract 



We present an efficient way to compute diagonal and off-diagonal n-point 
correlation functions for quantum spin-systems within the loop algorithm. We 
show that the general rules for the evaluation of these correlation functions 
take an especially simple form within the framework of directed loops. These 
rules state that contributing loops have to close coherently. As an application 
we evaluate the specific heat for the case of spin chains and ladders. 
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I. INTRODUCTION 



Numerical investigations of strongly correlated electron systems JTJ gained considerable 
importance in the last decade. The evaluation of non-diagonal correlation function and 
dynamical response function plays a major role in the context of correlated electron systems 
Jl|fJ]. On the other hand, there are only very few investigations of non-diagonal and/or 
higher-order correlation function in the context of quantum spin-systems. Indeed, it has been 
realized only recently, that non-diagonal correlation function might be calculated efficiently 
within the loop-algorithm |J . The loop-algorithm |4[] has established itself as the method of 
choice for quantum-Monte Carlo (MC) simulations of non-frustrated quantum spin systems. 

The key observation here is the fact, that local updating dynamics in a MC simulation 
creates strongly correlated configurations for gapless quantum spin systems at low temper- 
atures. Since the samples are then not statistically independent, the statistical error bars 
do decay only very slowly with the number of samples. One way to state this problem is to 
say, that the autocorrelation time r auto for the samples of spin-configurations created with 
the MC-walk increases (in generally exponentially) at low temperatures. 

Most efficient MC procedures implement consequently global update dynamics. Exam- 
ples of these procedures are the clusters algorithms ||. Designed to circumvent the critical 
slowing down, these methods have been intensively used to study classical statistical systems 
near critical points, where the problem of large r auto is very severe. 

The loop algorithm || can be considered as a generalization of classical clusters al- 
gorithms to quantum models. In fact, it gives a prescription on how global updates can 
be performed in quantum systems. As we will see this prescription lays on the geometric 
interpretation of the transformation from a quantum system to a statistical model of ori- 
ented loops. The MC procedure can be implemented then directly on the loops. It has the 
advantage that the updating dynamics defined on the loops generates statistically nearly 
independent configurations. The autocorrelation time is therefore about just MC step and 
the corresponding operators can be measured at every MC step avoiding both 'waiting times' 
and substantial increments of the variance (statistical error bars). 

In addition, a loop has another remarkable property; starting from an allowed spin con- 
figuration, constructing a loop and then flipping all spins in one loop (flipping the orientation 
of the loop) one obtains a new allowed configuration. This observation allows to compute 
the expectation value of operators not only in one configuration per MC step but in all 
configuration related to it by flipping any number of given loops. This procedure is usually 
called improved estimator [||]. 

The purpose of this work is to extend the algorithm to the computation of higher order 
(and non-diagonal) correlations functions. As we will see it involves dealing with two or more 
loop contributions. In particular we will focus on the specific heat cy, which, in the past, has 
been considered a major challenge for Monte-Carlo simulations |7|]. We will show, that the 
direct evaluation of the higher-order (non-diagonal) correlation functions contributing to Cy 
allows for improved estimators and such to gain one order of magnitude in computational 
efficiency. The method that we presented is valid in any dimension. 
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II. THE LOOP ALGORITHM 



A nice review of the loop algorithm can be found in Ref. ||. Here we start with a 
short introduction in order to introduce the notation used further on for the evaluation of 
higher-order correlation functions. 

The loop algorithm is most easily understood in the checkerboard picture for a discrete 
number of Trotter slices Nt] the generalization to continuous Trotter time || is straight- 
forward. This picture, which is based on the Suzuki- Trotter decomposition, describes in a 
graphical way how the interacting spin system wave function evolves in discrete imaginary 
time. 

The Suzuki- Trotter formula |J maps a quantum spin system in dimension d onto a 
classical spin in dimension d+1. The partition function of the original quantum spin model 
is hereby written in terms of the trace of a product of transfer matrices defined in the 
classical model. 

To illustrate the method we consider an inhomogeneous one-dimensional XXZ model 
H = Hi + H 2 on a bipartite chain of length L: 

i=2m j=2m+l 



jXY 

Hi = ^— {sfs i+1 + S i Sf + -^j + jf sf ' sf +1 , 

where the sign of the term ~ Jf Y has been choose to be negative by an appropriate rotation 
of the spins on one of the two sublattices. This is always possible on a bipartite lattice and 
allows for positive transfer matrix elements (absence of the sign problem). The decomposi- 
tion H = Hi + H 2 allows for the use of Totter-Suzuki formula M for the representation of 



the partition function Z = Tr exp(-^-H) 



N T ' 



N T 

Tr II 5>£2l exp(-ArlTi) ex P (-ArH 2 )\^} + 0(Ar 2 ) 



n=l «n 



where At = (3/Nt- Here we have introduced representations of the unity operator 
Xi Qn \(f>a ){4>a I m between any of the imaginary time slices. 

Since H\ and H 2 are sum of local operators that commute with each other, we may write 
the wave function as the product of the local basis in say z-component of spin, = ®i|cr?), 

with Oi =|, |. 

In the checkerboard lattice the interaction between two consecutive pairs of spins is 
graphically denoted by shaded plaquettes (see Fig. |1|). There are two spins interacting per 
plaquette so a 4x4 transfer matrix Tj can be defined in each plaquette, which depends only 
one the coupling constants. 
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FIGURES 



FIG. 1. Illustration of a plaquette with a spin flip process which corresponds to the matrix 
element (fj |exp(— AtH)\ J,|) of the transfer matrix. 
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The partition function Z is then, up to terms order 0(Ar 2 ), the trace of a product of transfer 
matrices: 

Z = Tr[exp(-/3#)]=Trn f ® ^ ( (g) T, 

n=l \i=2m / \j=2m+l 

As a next step beyond this standard representation of d- dimensional quantum models in 
terms of classical statistical systems [|T(J we expanded the transfer matrices T = TJ 7 p^ 'Af^ 

in terms of certain matrices such that the weig ht p{ y) > are non-negative. This is, 
in general, not possible for all models. For the XXZ with Jf y > Jf we can choose: 
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where p\ l} = ±(exp(-Ar Jf /2) + exp(-Ar Jf y /2)) exp(Ar Jf /4), pf = ±(exp(-ArJ?/2) - 
exp(-Ar y /2)) exp(Ar Jf /4) 

and pf ] = |(-exp(-ArJ?/2) + exp(ArJf y /2)) exp(ArJf /4). We then obtain for the 
partition function 



n=l i=2m \ 7 / j=2m+l V 7 / 



(1) 
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FIG. 2. Evolution of worldlines of up and down spins in imaginary time. Periodic boundary 
conditions are assumed in both space and imaginary time. Note that we define worldlines for both 
up- and down-spins, which differ by the direction in imaginary time. 



Eq. ([]]) can be interpreted in a geometrical way (see Fig. |2|). In the checkerboard pic- 
ture the M^ 1 ' matrices can be understood as different ways in which the worldlines can 
be broken in every plaquette and are usually called breakups. By taking one breakup per 
every plaquette we force the worldlines into closed paths which we call directed loops (see 
Fig. |3]). A directed loop therefore follows the worldline of an up-spin when it evolves in 
positive Trotter-time direction and the world-line of a down-spin when it evolves in negative 
Trotter-time direction. 

In Fig. |] we show the graphic representation of the breakups The lines now 

represent the directed loop segments. 




FIG. 3. Illustration of loop breakups, the directed lines represent the loop segments. From left 
to right the vertical (M^), diagonal (M^) and the horizontal breakup (M®) are shown. 
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Eq. ([!]) states that the partition function can be obtained as a sum over all breakups. 
As a sum over all breakups is equivalent to a sum over all loop configurations {1} we may 
rewrite Eq. ([]]) as 

Z = £>({/}) Tr I] (g) AfM (g) , (2) 

{2} n=l i=2m j=2m+l 

where p({Z}) = UiP^ TljP^- Eq. © leads to a very efficient MC-algorithm [|J: (a) Choose 

loop-breakups with probabilities p^ 1 ^. (b) Construct the loop configuration {1} and 

flip all loops with probability 1/2. (c) Measure any desired operator in all 2 Nl ^ 1 ^ spin 
configurations reachable with independent loop flips (improved estimators), where Nl({1}) 
is the number of loops in the loop configuration {I}. 

For later use we rewrite Eq. (0) in a form of traces over individual loops. Noting that 
vertical and diagonal loop segments do not change the spin-direction (see Fig. |3p, we may 

associate the 2x2 identity matrix o~° = ^ J ^ ^ with vertical and diagonal loop segments. 

As horizontal loop segments do change the spin-direction, we associate the Pauli-matrix 

1 

1 



° x = ( i n ) with them. We then may rewrite Eq. (0) as 



^=emw) n n»n* 7 " > ( 3 ) 

{;} ie{i} * 

where /x is an index running over loop I and 7 M = 0,x. Tr; denotes the trace over loop 
I. Since TiiY[ fl cr y ' 1 = 2, Eq. (H) is equivalent to a statistical mechanical model of oriented 
loops, Z = E W P({1}) 2^«'». 

III. CORRELATION FUNCTIONS, IMPROVED ESTIMATORS 

The expectation value of an operator (9 is 

(O) = Tr(e>exp(-/?F)) = ^(0 Q |C?|^)(^| exp(-/3^)|0 Q ) (4) 

If (9 is diagonal in the basis then this procedure is straightforward. The updating 

procedure generates a sequence of configurations q mc {iuc = 1 ■ ■ ■ Nmc), according with 
the distribution function of the system. In these configurations O takes a well defined value 
0(ci MC ), therefore: 

(0) = ^-T, 0{c iMC ) . (5) 

The loop algorithm allows to measure an operator not only in q mc but in all config- 
urations related by loop flippings. We illustrate the use of these improved estimators by 
computing O = S^S* (here indices x and y label both space and Trotter time (see Fig. 
H). When x and y belong to different loops the orientation can be changed independently 
and the total contribution cancels. By the contrary when x and y are on the same loop 
the orientations of the loop in both sites are linked and these terms contribute for the two 
possible orientations of the loop. 
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FIG. 4. Summing over all possible loop orientations, the two-loop contributions cancel each 
other for S^.Sy. Only when the two spin-operators act on the same loop we get a non- vanishing 
contribution. 



We will consider now the problem of non diagonal operators. The expectation value of 
a non diagonal operator O' in the loop picture is, see Eq. (0): 



{&) = ]T p{{l}) TrT M (7i) (g) 



(6) 



{0 



n=l i 



where T() means proper imaginary time ordering. Let us take as an example the two- 
point correlator O' = S^S~ , 





FIG. 5. The action of S~S+ can be represented as the insertion of a new kind of plaquette, 
here depicted in black, which acts on two loop segments. This operator flips the spins and changes 
therefore the orientation of the two loops for the remaining segments. Left-picture: The arrows 
denote the spin-direction. Right-picture: The arrows denote the direction of the loops. 
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Graphically the evaluation of an operator can be interpreted on the checkerboard frame- 
work as the insertion of a new kind of plaquette. In Fig. [5] we show the action of that operator 
in the checkerboard picture. We note that an off-diagonal operator in general reverses the 
direction of one or more loops. The loop configurations generated by the MC updating- 
procedure does, on the other hand, only generate loops with well defined loop orientations. 
Nevertheless there is a close connection between these two types of configurations which is 
easy to understand in graphical terms. 





FIG. 6. If S£ and S~ operate on different loops none of them can be closed consistently in 
terms of loop orientation, represented here by arrows. When both operators act on the same loop 
that configuration contributes to (S£S~). 



In Fig. |5] it is shown how the flipping of one spin 'propagates' through the loop, changing 
the orientation of the loop from that point. Thinking in terms of oriented loops it is obvious 
that with only one of these flipping processes (S£ or S~) per loop, it is not possible to 
close the loop consistently. To reestablish the original loop orientation it is necessary to 
have an even number of properly ordered S~ or S + operators on the same loop to close it 
consistently in terms of loop orientation variables. A loop which is not properly closed does 
not contribute to (O f ). Then we can establish that for a two-point correlation function we 
only obtain a contribution when x and y belong to the same loop (see Fig. ||). Eq. (f|) could 
suggest that measurements of non diagonal operators consume more computing time than 
diagonal operators, but using this graphical picture we note that both computations can be 
implemented in an equivalent way. 

These ideas can be justified in formal terms using Eq. (Q) and Eq. @. The and 
S y operators are placed in between of two a 7 matrices belonging to neighboring plaquettes 
and traces can be taken again independently in each loop. We define the 2x2 matrices 

a + = ^ |j q ^ and cr~ = ^ ^ °\ For positive loop-direction (with respect to the Trotter 

direction) S + is equivalent to a + , for a directed loop segment with negative loop-direction 
S + is equivalent to a~ . For S~ it is just the other way round. The loop direction of relevance 
here is the one before the insertion of either a S + or a S~ operator. 
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We start considering contributions to (S+S y ) where the loop-direction at site x is up 
and down at site y (see Fig. H). The expectation value of the non-diagonal operator S£S~ 
then becomes (compare Eq. (|3|)) 

(S+S y ) - Wp({l}) T lata; J] Tr, n^ 7 " ) • (7) 

« \ im » / 

Here T means proper time and space ordering. When cr+ and <jy are placed in different 
loops, the traces taken in these two loops cancel. If they are in the same loop the trace 
taken in that loop equals 1 (and not 2), independently of the spin-configuration. We will 
prove this last point now. We start by writing the partial trace of the loop containing <7+ 
and <7 y as 

where we neglected the cr° matrices, as they are just the identity matrices. We note that z\ + 
z 2 is even since {<J X ) 2 = cr° and because we are considering a loop which did contribute to the 
partition function Z before the S£S~ operators were inserted. The a x matrix corresponds 
to a horizontal loop segment and such to a change in loop direction. z\ needs therefore to 
be odd (and therefore also z 2 ), since one needs an odd number of directional inversions to 
arrive to a negative loop direction at site y, starting from a positive direction at site x. We 
may therefore rewrite T^ ++ ) (using again (<J X ) 2 = a ) as 

T (++) = Tr ; a+cr x cr+a x = 1 , 

as one can easily evaluate. Similarly one can consider the case when the initial loop directions 
are both positive at sites x and y. The expectation value of the non-diagonal operator S£S~ 
becomes then in this case 

(S+S y > - W P ({1}) T L$a~ I] Tr, n^ 7 " ] ■ (8) 

m \ im * J 

The corresponding one-loop contributions then have the form 

T< + -> = Tr z a+ (a x r ° y Kf 2 = Tr* a+a y = 1 , 

since both z\ and z 2 have to be even in this case. Similarly one can consider the two 
remaining cases of loop directions down/up and down/down at the sites x and y. It is 
worthwhile noting, that one easily proves along these lines the expected result (S+S^) = 0. 

IV. GENERAL CASE N-POINT CORRELATION FUNCTIONS 

In the last section we have shown how the loop orientation is the fundamental variable 
to deal with the computation of correlation functions using improved estimators. In fact the 
problem of n-point correlation functions can also be reduced to the study of how the loop 
orientation is changed by the action of some operators. 
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FIG. 7. Schematic example of two disconnected one- loop contributions to (O"). The dotted 
line in between two operators illustrates the case of two operators at the same Trotter time. 



We illustrate the case of two-loop terms for the four-point correlation function O" = 
S£ S~ S£,S~,. Here we consider the case relevant for the specific heat were (x, y) and (x', y') 
are pairs of real-space nearest neighbor (n.n.) sites at the same Trotter time. This oper- 
ator can generate several different kinds of contributions. The first one is the case of two 
disconnected one-loop contributions (see Fig. [?[). This is the case if 5+ and S~ act in one 
loop and S£, and S~, in a second loop. A second contribution arises if and S~, act in one 
loop and S~ and S£> in a second loop (see Fig. [|). We call this contribution a connected 
two-loop contribution. A third contribution arises when all four sites act on the same loop. 




FIG. 8. Schematic example of connected two-loop contribution. Two n.n. operators at the 
same Trotter time are connected with a dotted line. 

The evaluation of a single off-diagonal four-point operator O" does not pose a problem 
within the loop algorithm. For the case of interest, the specific heat a few additional points 
need to be kept in mind. The specific heat cy is given by 
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E<(s x -s y )(s x ,-s y ,)}- £<s x -s y > 



(9) 



where, again, (x, y) and (x', y') are pairs of (real-space) n.n. sites on the Trotter lattice. 
The first term of Eq. (^j) is a local energy-energy correlation function. When, x and y 
belong to a loop and x' and y' to another, we generate two-loop disconnected terms (as the 
one illustrated in Fig. ^) that can be computed from the expectation value of the internal 
energy, the second term of specific heat. The energy in a given MC-configuration, E iMC , can 
be written as a sum of the energy in the N L {i MC ) loops in this MC-configuration: 



E, 



IMC 



N L (i MC ) 

1=1 



With this definition we obtain 



(ind) 
Cy 



N 



MC IMC l^k 



k 

IMC 



N 



MC 



E 



J2 E Lc) -J2( E Lc) 



where cy = c v +c 



(conn) (ind) 



V 

the off-site terms, c 
on-site terms, Cy n \ 



(off) 



For the evaluation of the connected term Cy° nn ^ one has to evaluate 



. ^ v , where the pairs (x, y) and (x', y') are disjunct, separately from the 



where they are not disjunct: Cy° nn ^ 



>//) 



S on ) 



Cy"" + Cy'"' . By spin-algebra 
the on-site terms reduce to general two-point correlation functions. The (connected) off-site 
contributions fall in three categories, depending on the number S z operators involved (four, 
two or zero). The contributions with four S z operators have one and two loop contributions. 
A connected term with two S z operators has no two-loop contribution. Every correlation 
with two S z operators has the form S^S Z S£,S~,. If the indices x and y are not in the 
same loop the two S z operators act in different loops and their traces cancel for the reason 
explained in section III. The same reasoning is valid for x' and y' with the operators S + 
and S~ . Finally, terms with no S z operators can have two loop contributions (see Fig. [Bp 
and also one-loop contributions when the S + and S~ are properly ordered along the loop 
to close the loop coherently in terms of loop orientation. On the left of Fig. ^| we see that 
an arbitrary insertion of the operators S + and S~ can produce a conflict on the orientation 
of the loop. Technically, the value of the trace taken along the loop will depend on the 
structure of the correlator. This structure determines the order of the insertion of the a + 
and a~ matrices. For example the trace along the loop on the left of Fig. [| is: 

Tr(a x a x a x ,a y a x ap) = . 

For the loop on the right of Fig. || it is: 

Tr(a x aV y a+,crVy,) = 1 . 
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FIG. 9. On the left we show an example of non-contributing configuration to the specific heat. 
The loop orientation is ill defined and therefore this configuration does not contribute. On the 
right we see a contributing one-loop configuration. 



It is possible to evaluate certain off-diagonal operators O by an alternative method. The 
condition is, that the operator can be expressed by a sum of local operators which do involve 
the same pairs of sites (I, I') as the Hamilton-operator H = Y^a,i') Hij>: O = Yla,i') ®i,v- It is 
then possible to compute (O) by a reweighting method. The idea is to extend the plaquette 
of the checkerboard representation by new internal degrees of freedom, J2p \4 > p){4 > p\ ( see Fig- 
T~0|) . The reweighted matrix element of (O x ^ x >) is then 



Q( n ) 



Jx,x') = 



a n ,a n +iV~> ~ l /A n )l I A u \| 



(10) 



where x and x' denote combined space-time indices. For a given spin-configuration 
q mc = {(j)^a}\( n = 1, ■ ■ ■ , Nt)} the off-diagonal expectation value of C(q mc ) is 0(q mc ) = 
I/NtE^Uu) 0S, an+1 (x,x') and (O) = l/N MC ^ MC 0(c lMC ) (see Eq. (§)). 



(n+l) 

a. 



4>, 




(n) 
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FIG. 10. New structure of the plaquette in the reweighting method. The grey plaquette is 
the conventional plaquette where the evolution in imaginary time takes place, the black plaquette 
represents the operator O on the basis of a z . The product of the two matrices generates a new 
composite plaquette where the new weight is defined. 

The reweighting method may also be applied to specific heat, which is the sum of products 
of local operators. 

From the point of view of the complexity of the algorithm, measuring four-point cor- 
relation functions requires more computing time than two-point correlation functions. For 
the latter is only necessary to know whether or not two sites are in the same loop. This 
information can be obtained at the same time the loop is constructed and consequently the 
computing time remains proportional to LNt- For n-point correlation functions the situa- 
tion is more complex. In this case, there are contributions involving two or more loops and 
at the same time non-diagonal operators give different contributions depending on how they 
are ordered on the loop. In practice this depends on the shape of the loops in each configu- 
ration. A rigorous study of the performance of the method must include an analysis of the 
behavior of the statistical errors as a function of the temperature, size, number and type 
of operator involved in the correlation functions and the details of the Hamiltonian. This 
detailed analysis of technical aspects of n-point correlations will be presented elsewhere. 

V. RESULTS 

As an application of the rules explained in this paper we have computed the specific heat 
for a Heisenberg chain and for a ladder with J± = 0.5J (which corresponds to the ratio for 
the ladder-compound Sr 14 Cu240 4 i |H§)- 
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FIG. 11. Specific heat for the Heisenberg model in the 8-site chain with J = 2, as a function of 
temperature in units of J. The two sets of data correspond to the QMC simulations with improved 
estimators (open squares) and with the reweighting method (filled triangles), for the same number 
of QMC steps. The solid line is the exact diagonalization data. 



In Fig. |TT| we compare exact diagonalization results with the results using the method 
described above and the reweighting method for the same number of MC steps. The error 
bars in these two methods are also compared. For the lowest temperature the error bar with 
improved estimators are 6 times smaller. Taking into account that error bars decay as ^ 
we expect that without using improved estimators 36 times more MC steps are necessary to 
get equal size error bars. The statistical errors are amplified by the factor (3 2 . This factor 
and the substraction of similarly large numbers lead to large error bars at low temperatures. 

In the Fig. O we present results for the specific heat of a 100-site Heisenberg chain. 
To reproduce the linear regime at low temperatures it is necessary to perform a careful 
extrapolation to At — > taking half a million of MC steps for each At values and 10 
different Nj< values ranging from 20 to 200. 
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FIG. 12. Specific heat for a 100 sites Heisenberg with J=2.0 chain using improved estimators, 
as a function of temperature, in units of J. The error-bars are smaller than the symbol sizes. The 



solid line is the exact Bethe-Ansatz result for the infinite-chain |11]. 



In Fig. [13] we present results for the two-leg ladder of 2 x 201 sites with twisted boundary 
conditions (i.e. for J± = this system corresponds to a L=402-site Heisenberg chain). 
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FIG. 13. Specific heat for the 2 x 201 ladder with twisted boundary conditions as a function of 
temperature in units of J. The values of the couplings are J = 1.0 and J± = 0.5J. 



VI. CONCLUSIONS 

We have presented detailed rules on how to evaluate general, off-diagonal n-point Greens 
functions within the loop algorithm. These rules have a very simple interpretation in the 
picture of oriented loops. They state that the loop-orientation has to close coherently when- 
ever a certain number of non-diagonal operators are inserted. We have shown how to apply 
these rules to the case of the specific heat and presented results for the lD-Heisenberg model 
and a ladder system. 
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